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Large  Space  Structure 
Charging  During  Eclipse  Passage 


1.  INTRODUCTION 

The  large  variations  in  spacecraft  potential  measured  relative  to  the  plasma 

that  have  been  observed1  are  believed  to  be  the  result  of  the  dominance  of  the  elec- 

2 

tron  flux.  This  may  be  explained  as  follows.  According  to  Garrett,  the  two 
plasma  populations,  protons  and  electrons,  can  each  be  characterized  as  a  pair  of 
Maxwellian  populations: 


each  of  which  has  a  number  density  (N)  and  a  Maxwellian  temperature  (T).  The 
energy  of  each  particle  is  proportional  to  the  mass  of  the  particle  and  the  square 
of  its  velocity;  thus,  for  a  given  energy,  protons,  being  three  orders  of  magnitude 
more  massive  than  electrons,  have  velocities  43  times  lower.  Temperature  and 
number  flux,  being  comparable  on  the  macroscopic  scale  to  energy  and  velocity  on 

(Received  for  publication  15  January  1980) 

1.  DeForest,  S.  E.  (1972)  Spacecraft  charging  at  synchronous  orbit,  J.  Geophys. 

Res.  77(No.  4):651.  - 

2.  Garrett,  H.  B.  (1977)  Modeling  of  the  Geosynchronous  Plasma  Environment, 

A FGL-TR -  77-0288,  AD  A053  164. 
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the  microscopic,  have  a  similar  relationship,  leading  to  the  dominance  of  the 
electron  flux. 

Under  the  influence  of  only  these  hot  Maxwellian  plasmas  (kT  greater  than 

2 

10  eV),  any  object  in  the  plasma  absorbing  charged  particles  would  collect  a  net 
negative  charge,  attaining  an  electric  potential  large  enough  to  repel  the  incident 
electron  flux  (-10  V  to  -20,  000  V).  However,  other  effects  can  change  the  situa¬ 
tion.  If  an  additional  dense  "cold"  plasma  population  (kT  less  than  10  eV, 

3 

N  greater  than  10  n/cm  )  is  present,  then  any  negative  voltage  is  neutralized  by 

3 

an  enhanced  ion  flux.  This  effect  suppresses  charging  in  low  earth  orbit  but  not 

at  geosynchronous  altitudes,  where  the  cold  population  is  extremely  thin 
-2  3 

(N  =  10  n/cm  ).  When  the  structure  is  illuminated  (for  example,  a  satellite  in 
the  sunlight),  electrons  are  removed  from  the  surface  by  the  photoelectric  effect. 
Often  this  photoelectron  current  dominates  the  ambient  electron  current,  resulting 

4 

in  a  small  positive  potential  (0.  1  to  10.  0  V). 

As  a  small  satellite  passes  from  eclipse,  its  voltage  will  undergo  a  drastic 
change  (~  1000  V  in  1  min),  corresponding  to  the  change  in  illumination.  If  it  is 
large,  as  mentioned  earlier,  the  gradient  of  illumination  will  lead  to  a  voltage 

differential  across  the  surface.  As  will  be  seen,  the  magnitude  of  this  differential 

-3  3 

potential  will  range  from  10  to  10  V  for  a  10-km  structure. 


2.  BASIC  THEORY 


If  Jp  is  the  total  current  density  to  a  small  spacecraft,  then1 


J  =  J  -  (J.  +  J  +  J  .  +  J,  )  +  J..  -  J  -j-2- 
p  e  i  se  si  bs  th  pe  Ag 


(2) 


where 

Je  =  the  ambient  electron  current  density 
J.  =  the  ambient  ion  current  density 

Jse  =  the  secondary  electron  current  density  caused  by  primary  ambient 
electrons 

Jgi  =  the  secondary  electron  current  density  caused  by  primary  ambient  ions 


3.  Whipple,  E.C.  (1965)  The  Equilibrium  Electric  Potential  of  a  Body  in  the 

Upper  Atmosphere.  NASA  X-615-65-296. 

4.  Garrett,  H.B.,  Pavel,  A.L.,  and  Hardy,  D.  A.  (1977)  Rapid  Variations  in 

Spacecraft  Potential.  AFGL-TR -77-0132,  ADA046  35UT 

5.  Rosen,  A.  (1975)  Spacecraft  Charging:  Environment  Induced  Anomalies. 

Paper  75-91,  AlAA  13th  Aerospace  Sciences  Meeting.  ’ 
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Jj_)S  =  the  backscatter  electron  current  density 

J  ^  =  the  "thruster"  current  density,  caused  by  ion  drive,  plasma  beams,  etc. 

Jpe  =  the  photoelectron  current  density 

A  =  the  sunlit  cross-section  area  of  the  spacecraft,  and 
c 

Ag  -  the  surface  area  of  the  spacecraft 

The  corresponding  values  for  Je,  Jge,  Jgi,  Jpe,  and  J^g  as  functions  of  volt¬ 
age,  number  density,  and  temperature  have  been  derived  by  Tsipouras  and 

g 

Garrett  (assuming  a  thick  plasma  sheath,  an  isotropic,  two-Maxwellian  plasma, 

and  an  aluminum  surface);  these  values  are  listed  in  the  Appendix.  The  thruster 

current  density  has  been  assumed  to  be  zero  for  this  study,  and  any  capacitance, 

inductance,  or  local  differential  charging  effects  are  ignored.  For  a  system  in 

£ 

equilibrium,  the  plasma  current  I  (that  is,  the  plasma  current  density  multi¬ 
plied  by  the  surface  area  A  )  will  be  equal  to  zero,  defining  by  Eq.  (2)  a  unique 
voltage  for  a  given  situation  that  will  be  the  potential  on  the  spacecraft  as  meas¬ 
ured  relative  to  the  plasma.  While  the  currents  are  time  varying,  the  time  con- 

7 

stant  for  Eq.  (2)  to  be  satisfied  is  on  the  order  of  milliseconds,  and  the  assump¬ 
tion  of  equilibrium  is  justified. 

The  extension  to  a  large  structure  such  as  a  space  based  radar  (SBR)  or  solar 
power  satellite  (SPS)  is  made  by  two  changes  in  Eq.  (2).  First,  the  ambient  plasma 
equations  change,  attenuating  the  number  flux  of  the  attracted  species  (thin  sheath 
approximation);  that  is,  the  current  of  the  attracted  species  is  assumed  to  be  inde¬ 
pendent  of  the  voltage;  there  is  no  change  in  the  repelled  species  (see  Appendix). 
Second,  an  ohmic  current,  the  current  gained  or  lost  to  other  parts  of  the  space¬ 
craft,  is  added  to  Eq.  (2).  Equation  (3)  now  gives  the  net  current  to  a  particular 
surface,  or  node  of  the  spacecraft: 


^n  V*  +  *ohm 


(Please  note  that  any  future  reference  to  Ip  includes  the  ambient  plasma,  second¬ 
ary  electron,  backscattered  electron,  and  photoelectron  currents;  similar  con¬ 
ventions  hold  for  Jp,  and  the  phrases  "plasma  current"  and  "plasma  current 
density."  Also,  direction  of  current  flow  is  taken  to  be  the  direction  of  electron 
flow,  and  potentials  are  measured  relative  to  the  plasma. ) 

6.  Tsipouras,  P. ,  and  Garrett,  H.  B.  (1979)  Spacecraft  Charging  Model  —  Two 

Maxwellian  Approximation,  AFGL-TR -79—0 1 53,  AD  A 077  907. 

7.  Rothwell,  P.  L. ,  Rubin,  A.G.,  and  Yates,  G.K.  (1977)  A  simulation  model  of 

time-dependent  plasma  spacecraft  interactions,  Proc.  of  the  Spacecraft 

Charging  Conference,  AFGL-TR -77-0051/NASA  TMX-73537. 
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-  the  total  current  to  the  node 

-  the  ohmic  current  to  the  node 

-  the  plasma  current  to  the  node  (Eq.  (2)) 


Since  the  structure  is  assumed  to  be  in  equilibrium  as  before,  the  current  to  each 
node  must  be  equal  to  zero.  Finding  the  set  of  voltages  for  which  this  is  true 
constitutes  the  solution  to  the  charging  problem. 


3.  RESULTS 

A  set  of  equations  defined  by  Eqs.  (2)  and  (3)  was  solved  by  program  SHADOW 
using  iterative  means  on  the  CDC  6600  system  at  AFGL,  Hanscom  AFB.  Figure  1 
presents  the  results  for  Eq.  (2)  (that  is,  a  single  node).  The  horizontal  axis 
represents  actual  potentials  observed  by  the  NASA  geosynchronous  satellite  ATS-5 
during  21  eclipses  from  1969  to  1970;  the  vertical  axis  represents  the  potentials 
predicted  by  SHADOW.  The  points  plotted  fall  close  to  the  identity  line,  and  have 
a  standard  deviation  of  about  1200  volts,  confirming  the  validity  of  SHADOW  in 
calculating  satellite  potentials. 


Figure  1.  Observed  ATS-5  Potentials  in  Eclipse  versus  Potentials  Pre¬ 
dicted  by  Program  SHADOW 
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Figure  2  is  a  plot  of  the  single  point  satellite  potentials  V  normalized  to  the 

8  ® 

shadow  potential  V  against  position  in  orbit  Xm>  This  parameter  Xm  is 

defined  as  the  minimum  altitude  of  the  central  ray  of  the  disk  of  the  sun  as  seen 

from  the  satellite  (Figure  3a).  The  model  used  to  describe  the  illumination  of  the 

spacecraft  as  a  function  of  X  ,  F(X  ),  assumes  no  attenuation  due  to  the  atmos- 

m  m 

phere,  an  earth  radius  of  6378  km.  and  an  apparent  solar  radius  of  185.  2  km  (see 
Appendix  for  formula.  Figure  2a  is  a  plot  of  the  observed  potentials  (normalized 
by  the  eclipse  potentials  V  )  as  the  satellite  ATS-5  entered  or  left  eclipse,  whereas 

Figure  2b  is  a  plot  of  the  potentials  predicted  by  the  program  SHADOW,  assuming 

2  “62 

a  photoelectric  current  density  of  0.4  nA/cm  (or  4  '  10  A/m  ),  the  plasma  con¬ 
ditions  observed  by  ATS-5,  a  roughly  spherical  probe.  The  shape  of  the  satellite 
is  important  in  determining  the  ratio  Ac/Ag  of  Eq.  (2). 


Figure  2a.  Potential  During  Eclipse  Pas¬ 
sage,  VD,  Normalized  by  the  Potential  in 
Eclipse,  Vmax,  as  a  Function  of  the  Min¬ 
imum  Ray  Path  Height,  Xm,  to  the  Center 
of  the  Sun.  Twenty-one  examples  of 
actual  ATS-5  and  ATS-6  plasma  observa¬ 
tions  during  eclipse  passage  are  plotted 


8.  Garrett,  H.  B. ,  and  Rubin,  A.G.  ( 1978)  Spacecraft  charging  at  geosynchronous 
orbit-generalized  solution  for  eclipse  passage,  Geophys.  Res.  Lettrs. 

5(No.  10):865.  - 
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Figure  2b.  Potential  During 
Eclipse  Passage  for  Predicted 
vo/vmax  Values  Using 
SHADOW 


Next  a  series  of  runs  was  done  Cor  several  SPS  (Solar  Power  Satellite!  and  SBR 
(Space  Based  Radar)  models.  Each  structure  was  broken  into  a  series  of  segments 
or  nodes,  and  the  potential  on  each  node  found.  The  structures  were  assumed  to  be 
10  km  across,  with  a  thickness  of  1  cm.  Each  structure  was  divided  into  a  series 
of  10  to  100  nodes,  each  with  the  same  width  and  approximately  uniform  illumina¬ 
tion  (Figure  3b).  The  structures  were  also  assumed  to  be  homogeneous  with  the 
surface  properties  (photoelectron,  secondary  emission,  etc. )  of  ATS-5.  The 
conductivity,  number  of  nodes,  and  plasma  parameters  were  varied  in  order  to 
test  the  numerical  sensitivity  of  the  program  (10  nodes  gave  results  accurate  to 
~2  percent  of  the  100-node  model).  A  wide  variety  of  information  was  extracted 
from  each  run.  To  list  a  few:  the  potential,  ohmic  energy  dissipation,  the  ohmic 
current  at  each  point  in  the  structure,  and  the  total  voltage  drop  and  energy  dissi¬ 
pation  across  the  structure. 

The  potentials  in  the  penumbra,  based  on  a  thin  sheath  approximation,  tend 
to  be  twice  the  magnitude  of  those  potentials  predicted  using  the  thick  sheath 
approximation  (Figure  4).  This  is  a  direct  result  of  the  attenuation  of  the  attracted 
species  in  the  thin  sheath  approximation.  Also,  the  potential  gradients  were 
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Figure  3a.  Illustration  of  Parameter  X  Measurement.  Note  that 
Xm  is  the  height  of  the  minimum  ray  path  to  the  center  of  sun  (S) 
above  the  center  of  earth  (E)  as  seen  by  the  satellite  (SAT) 


Xj  (km) 

-5  0  +5 
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ELECTRICAL  ANALOG 

f  igure  3b.  Solar  Power  Satellite  Schematic  Showing  Coordinate  System 
and  Electrical  Analog  for  Simulation  Purposes 
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VMAX  f0*  P0INT 


Figure  4.  Plot  of  Eclipse  Potentials  Determined 
for  21  ATS-5  Events  Employing  a  Planar  (thin 
sheath)  or  Point  (thick  sheath)  Approximation 

greatly  affected  by  the  ohmic  currents  for  low  resistance  models,  differing  by  as 
much  as  500  V  from  the  high  resistance  potentials.  The  potentials  for  a  poorly 
conducting  surface  are  more  positive  toward  the  high  Xm  end  (sunlit)  than  for  a 
good  conductor,  whereas  those  toward  the  low  Xm  end  (shadowed)  tend  to  be  more 
negative.  This,  of  course,  is  due  to  the  characteristic  uniform  potential  of  a  per¬ 
fect  conductor.  This  latter  potential,  measured  relative  to  the  plasma,  can  be 
easily  determined  by  using  the  average  value  of  the  photoelectron  current  density 
J  •  F(X  )  across  the  structure  in  a  single  point,  thin  sheath  potential  calculation. 

P  The  results  of  calculations  using  ATS-5  data  from  Day  289,  1969  (a  particu¬ 
larly  severe  day  for  charging)  at  XQ  =  -150  km  are  plotted  in  Figures  5  to  10. 
Figure  5  shows  the  potential  drop  relative  to  the  shadowed  end  of  the  SPS  for  five 
different  resistances  (see  Table  1  for  plasma  parameters).  In  Figure  6  these 
potentials  have  been  normalized  by  dividing  by  the  total  drop  between  the  ends  of 
the  structure  AV.  The  high-resistance  curve  is  nearly  linear,  owing  to  the  fact 
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Figure  5.  Potential  Drop  Across 
a  10-km  Solar  Power  Satellite  with 
Its  Center  at  XQ  =  150  km.  Vari¬ 
ous  values  of  the  conductance  have 
been  assumed  while  the  ambient 
conditions  were  those  measured  by 
ATS-5  on  Day  289,  1969.  The 
dashed  line  is  a  similar  plot  for 
the  Space  Based  Radar 


Figure  6.  The  Potential  Di  ops  in 
Figure  5  are  Redrawn,  Normalized 
by  Division  by  the  Total  Potential 
Drop  across  the  Solar  Power  Satel¬ 
lite  to  Show  Variation  in  potential 
Across  the  Structure  Mot  e  Clearly 


Table  1.  Space  Plasma  Parameters  From  ATS-5 


Ent  ity 

Day  289, 
N  (cm-3) 

1969 

% 

T  (eV) 

Day  291, 
N  (cm  3) 

1970 

T  (eV) 

Electrons 

0.292 

1060 

0.  835 

998 

0.  550 

6070 

0.  227 

3910 

Ions 

0.057 

103 

0.  276 

319 

0.  690 

8090 

0.  911 

9530 
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Figure  7.  Maximum  Potential  Drop  Across  a 
Solar  Power  Satellite  Located  at  150  km 


that  F(X  )  is  linear  to  the  first  order  over  the  length  of  the  structure  at  ~  150  km. 
The  low-resistance  curves  show  a  definite  S-shape.  Also  plotted  in  Figure  5  are 
potentials  for  a  circular  10 -km  diam  structure  (SBR  model)  compared  to  the 
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square  SPS.  This  curve  has  been  compared  for  a  resistivity  of  10'  n  -m  for  the 
same  conditions  as  the  SPS.  There  is  little  difference  in  this  and  the  correspond¬ 
ing  SPS  curve;  however,  the  potential  drop  of  the  SBR  is  less,  most  likely  a  result 
of  the  reduced  SBR  area  compared  to  the  SPS. 

The  effect  of  resistance  on  the  total  voltage  drop  is  plotted  in  Figure  7.  The 
vertical  axis  is  the  voltage  drop;  the  horizontal  axis  is  the  common  logarithm  of 
resistance  in  ohms  for  a  10-km  structure  at  Xm  -  -150  km,  again  for  Day  289  of 

1969.  The  change  from  a  conductor  to  nonconductor  occurs  for  resistivities  of 
0  5 

10  to  10  D -m,  roughly  the  resistivity  of  semiconductors. 

7 

Figure  8  is  a  plot  of  the  voltage  drop  against  X  for  a  resistance  of  10  D  -m 
for  the  plasma  conditions  on  Day  289,  1969  and  Day  291,  1970,  and  assumes  a  point 
or  plane.  Figure  9  is  the  potential  drop  against  the  maximum  potential  in 
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Figure  8.  Maximum  Voltage  Drop  Across 
the  Solar  Power  Satellite  (and  Space  Based 
Radar)  for  Various  Plasma  Conditions  and 
Sheath  Assumptions 
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Figure  9.  Maximum  Voltage  Drops  at 
XQ  =  -150  km  for  the  21  ATS-5  Plasma 
Events  as  a  Function  of  V  ,  the 
Eclipse  Potential 
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L06(0[yO(°HMs-m)] 

Figure  10.  Ohmic  Power  Loss  in  Solar  Power  Satel¬ 
lite  as  a  Function  of  Resistivity.  The  Solar  Power 
Satellite  is  at  X0  =  150  km;  plasma  conditions  are 
for  Day  289,  1969 

the  earth's  shadow,  for  the  21  ATS-5  plasma  conditions  and  for  Xm  =  -150.  The 
points  are  contained  by  the  lines  Vprop  =  0.  06  V  +  120  and  V^rop  =  0.  03  V. 

As  Vmax  approaches  zero,  the  electron  current  is  dominated  by  the  photo¬ 
electron  current.  Thus  it  appears  that  at  V, _ =  0  V  there  would  be  a  minimum 

voltage  drop  dictated  by  the  photoelectron  current;  that  is,  the  differential  photo¬ 
electron  currents  by  themselves  would  drive  ohmic  currents  which  would  deter¬ 
mine  the  voltage  drop.  In  practice,  because  V  is  zero,  one  end  must  have 
V  >  0  and  the  photoelectron  current  becomes  negligible.  This  explains  the  region 
of  Figure  9  where  the  potential  —  0  but  the  voltage  drop  is  >100  V.  The  last  plot. 
Figure  10,  is  of  the  ohmic  power  dissipated  through  a  structure  as  a  function  of 
the  logarithm  of  resistivity.  The  value  of  resistivity  for  which  the  power  is  max¬ 
imum  has  been  empirically  fit  to  a  function  of  length  and  photoelectron  current: 
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Jp  =  photoelectron  Flux  (A /m  ) 
p  =  ohm-m 


4.  THEORETICAL  FIT 

Several  of  the  preceding  results  can  be  theoretically  explained  as  follows. 
Assume  that  a  homogeneous  plane  structure  of  uniform  thickness  T  and  resistivity 
p  is  divided  along  a  line  of  equal  illumination  into  two  sections  (Figure  3b).  The 
current  through  the  interface  is  equal  to  the  plasma  current  at  equilibrium  to 
either  section 


/vv«x)dA  =  as  (4) 

Here  J  (V  X)  dA  is  the  differential  current  from  the  ambient  environment  to  sur- 
’  P 

face  dA  at  potential  V  for  J  (V  )  =  0  =  J  (V  )  +  J  (V  );  X  is  a  coordinate  system 
X  1  x  p  x  o  x 

measured  perpendicular  to  the  interface  and  relative  to  the  center  of  the  structure, 
increasing  in  the  direction  of  increasing  illumination  (Figure  3b);  R  is  the  resis¬ 
tance  coordinate,  increasing  in  the  direction  of  increasing  X  (defined  as 
R  =  /p /A  dA).  Since 

dA  -  L(X)  dX  (5) 

p  dX 

dR  =  TOT  (6) 

where 

L(X)  is  the  width  of  the  structure  at  X, 


r 

J  jt(v,  x)  l(x)  dx  =  ^  1LM2D 

-S/2 


(7) 


S  is  the  length  of  the  structure  parallel  to  the  X  axis;  Xj  is  the  value  of  X  at 
the  interface. 


17 


Differentiating, 


C 


JT(V,  X)  1(X)  .  fa  g) 


This  may  be  rewritten  as 


-fjr  JT(V,  X)  =  D2V  +  g(X)  Dx  V 
where 


(9) 


g(x) 


=  dL(x)  1 
dx  L(x) 


(10) 


The  left-hand  side  of  Eq.  (9)  may  be  approximated  by  a  Taylor's  series: 


JT(v,  x)  •=  JT(vlt  x)  +  j^Vj.  xxv  -  vx)  +  Jy(vv  x)  — 


(11) 


Because  of  the  relative  insignificance  of  the  higher  order  derivatives  of  JT  (calcu¬ 
lated  numerically): 


JT(V,X)  i  10~10  A/Vm2 


JT(V,X)  1  10'14  A/V2m2 


JT(V,X)  1  10'18  A /V3m2 
dVJ  1 


and  so  on, 


all  but  the  first  two  terms  may  be  ignored  for  a  result  accurate  to  within  firs' 
order  for  |v  -  Vj  |  <  103  V 


JT(V,X)  =  JT(Vj,X) 


+  Hv  JT(V1-X)  <V  -  V1> 


(12) 


(13) 


Substituting  in  the  value  of  V  for  which  is  zero,  namely  V^,  we  find: 

JT<V'X)  1  dV  Jt(V2*X>  (V  ‘  v2> 

Similarly,  as  Vz  is  a  function  of  X  and 

h  Vz(X)  S  10-1  V/m 
2 

~2  VZ(X)  '=  10'3  V/m2 
ctX 

and  so  on. 


VZ(X)  =  Vz(0)  +  Vz(0)  X 

Thus 

JT(V,X)  =  Jx(v2.  X)  (V  -  Vz(0)  -  Vz(0)  X) 

»» 

Because  J-j-(Vz>  X)  is  insignificant  for  the  range  considered, 
J^,(Vz,X)  ’=  JT(Vz(0),  0)  , 

jt(v,  x)  =  JT(vc,  0)  (V  -  vc  -  Vz(0)  X) 
where 

Vc  =  Vz(0)  ,  center  of  structure. 

Substituting  Eq.  (17)  into  Eq.  (10), 

T  JT(Vc-  0)  (V  •  Vc  ’  Vz(0)  x>  =  dx2V  +  (DXV)  g  (x) 
or,  putting  all  references  to  V  on  the  same  side  of  the  equation, 

-  ^  JT(Vc,  0)  <VC  +  Vz(0)  X)  =  (Dx2  +  g(x)  Dx  -  J^.(Vc,  0))  V 


(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 
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n 


we  find  that  Eq.  (20)  can  be  solved  using  standard  numerical  techniques.  However, 
if  the  width  of  the  structure,  L(X),  is  a  non-zero  constant  (SPS  rectangular  model), 
then  g(X)  is  zero,  and  Eq.  (20)  can  be  algebraically  solved  for  V  as  a  function  of 
X: 

c.  (a-X)  Cn  (— a«X) 

V  =  Te  +Te  t-agX  +  aj  (21) 

where 

Vc 

V>) 

<j;(Vc,  OJp/T)1/2 
Assuming  that  V(0)  =  V 

V  =  c%  sinh  (a3X)  +  a2X  +  &x  (22) 

By  Eq.  (7)  when  X  =  -  s/2,  dV/dX  =  0,  and  by  Eq.  (22), 

g  =  Cj  a3  cosh(a3X)  +  a2  , 
so  substituting  -  s/2  for  X,  we  get  a  value  of 


_ a2 

C1  "  a3  cosh  (-sa3/2) 

Thus,  the  final  equation  for  the  SPS  model  is 


(23) 


V  =  a^  sinh  (a3X)  +  a2S  + 


where 


a 


1 


V 


c 


V'z(0) 

(J^.(Vc,0)p/T)1/2 
-  a2/a3  cosh  (-sa3/2) 


(24) 


This  equation  fits  the  results  of  the  numerical  calculations,  program 
SHADOW,  to  within  approximately  10  percent.  However,  to  achieve  this  accuracy, 
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it  is  necessary  to  have  values  of  and  for  the  plasma  data  used,  and  to  stay 
within  the  ranges  specified  by  the  approximations  used: 

s  <  XL 

where  XL  is  the  lowest  value  of  XM  for  any  point  on.the  structure,  and 

V  (X)  <  0  for  all  X  within  the  structure  . 
z 

The  worst  instances  of  charging  were  observed  using  data  from  Day  289,  1969, 
and  from  these  data,  values  of  10  *  V/m  and  10  A/Vm  for  and  JJj,  were 
derived.  These  values  may  be  taken  to  be  a  worst  case. 


5.  CONCLUSION 

The  results  of  this  study  indicate  that  potential  differences  of  1,  000  V  across 
10-km  structures  are  possible  under  even  the  mild  plasma  conditions  and  with  the 
relatively  weak  photoelectric  current  of  the  ATS-5  eclipses.  Furthermore, 
engineers  must  be  made  aware  of  the  effects  of  surface  properties  and  resistivities 
of  various  spacecraft  materials  on  spacecraft  charging.  One  of  these  effects,  pre¬ 
viously  unsuspected,  is  the  existence  of  a  resistivity  for  which  the  energy  dissi¬ 
pation  through  a  structure  becomes  a  maximum  and  which  is  a  function  of  length 
and  photoelectron  current  density. 

Four  tools  are  provided  the  engineer  for  calculation  of  the  voltages  across  a 
large  space  structure  passing  into  eclipse:  (1)  the  program  SHADOW,  which  uses 
a  detailed  but  efficient  numerical  algorithm;  (2)  differential  Eq.  (9)  which  may  be 
solved  using  numerical  methods;  (3)  Eq.  (20),  which  can  be  solved  with  greater 
efficiency  but  with  less  accuracy;  and  (4)  specifically  for  the  SPS  model,  Eq.  (24), 
which  provides  results  with  extreme  efficiency  and,  so  far  as  the  limitation  on 
shape  takes  it,  the  same  accuracy  as  Eq.  (20). 

Because  these  methods  ignore  the  effects  of  capacitance  and  inductance,  they 
must  be  considered  qualitative.  However,  knowledge  of  the  potential  and  thus 
charge  distribution  can  be  used  to  estimate  the  former  effects.  Therefore,  these 
techniques  may  prove  an  invaluable  aid  in  predicting  the  charging  effects  on  large 
space  structures  passing  into  eclipse. 
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Appendix  A 

Current  Values  as  Functions  of  Voltage, 
Number  Density,  and  Temperature 


CURRENT  DENSITY  EQUATIONS 

For  incident  ambient  electrons  (or  ions)  the  current  density  is  the  sum  of 
the  current  densities  of  two  populations  where  the  current  density  for  each 
population  is: 

J  =  JQ  exp  ( - 1  q V |  /kT)  (for  the  repelled  species) 

J  -  JQ  (1  +  |qV|  /kT)  (for  the  attracted  species) 

where  is  the  population's  current  density  at  V  =  0,  q  the  charge  of  the  particular 
species,  and  kT  the  Maxwellian  temperature  of  a  single  population.  The  current 
densities  of  the  secondary  electron  and  backscatter  electron  fluxes,  known  collec¬ 
tively  as  the  secondary  current  densities,  are  the  sum  of  the  secondary  current 
densities  caused  by  each  of  the  primary  populations;  that  is,  the  secondary  elec¬ 
tron  current  density  is  the  sum  of  the  secondary  current  densities  caused  by  each 
of  the  two  electron  populations  and  the  two  ion  populations,  and  the  backscatter 
current  density  is  the  sum  of  the  backscatter  densities  of  the  two  ambient  electron 
populations.  On  this  population  by  population  basis,  the  equations  are  relatively 
simple: 


25 


Jsel=Jei*a(V-Tel> 


Jsil  =  Jil*b<V’Til) 


Jbsl  -  Jel*c(V'Tel> 


and  so  on. 


where 


a(V,  Te)  =  0 

V  >  0 

=  0.  8431  +  (-0.  8424)  exp  (-0.  0286*Te) 

V  <  0, 

T  <  200  eV 
e 

a(V,  T  )  =  0.  1431  +  0.  665  exp  (-0.  0003*Tg) 

V  <  0, 

T  >  200  eV 
e 

c(  V,  T  )  =  0.  222 

V  <  0 

=  (0.  263  +  (-0.  0218)  exp  (-0.  00901*V)) 

+  (0.  11536  +  (-0.  00328)  ■  exp  (-0.  15127*V)) 

•  exp  ((-0.  00081  +  0.  00019  exp  (0.  02217*V))*Te)  V  >  0 

b(V,  T()  =  0  V  >  0 

b(V,Ti)  =  (4.  78  -  0.  6526  exp  (V*0.  9783*10~4)) 

+  (-0.  5299  -  3.  78  exp  (V*0.  1844*10'3)) 

•  exp  (T. *( -0.  5715*10-4  -  0.  1829*10‘3 

•  exp  (0.  1447*10'3*V)))  V  <  0 

These  values  of  a,  b,  and  c  have  been  calculated  for  aluminum.  They  have  been 
multiplied  by  arbitrary  coefficients  to  fit  the  ATS-5  observations  (that  is, 

Jse  V  a-  X  =  1'3:  Jsi  =Y  '  V  b-  Y  =°-55:  Jbs  =  Z‘  Je‘  c>  Z=  °-4)‘ 
The  photoelectron  current  density  is: 


=  J  = 

*F(X  ) 

V  <  0 

peo 

m 

=  J 

peo 

*F(X  )*(1  +  V/0.  7)"2 
m  ' 

V  >  0 
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when 


J  is  the  saturation  photoelectron  current  (V  =  0,  X  -  200) 

peo  m 

F(X  )  is  the  fractional  illumination  of  the  spacecraft: 

2 

^  (A  -  sin  A)  +  (B  -  sin  B) 

m'  ~  1  ~ 


where 


A  =  2  cos 


(<Xm  +  RE>2  +  RE-R^ 

V  -2(Xm  +  RE)RE  -  / 


B  =  2  cos 


/  (Xm  +  RE)2  +  Rs-RE 

\  2(X  +  R_)  R 

\  m  E  s 


Re  =  earth  radius 

Rg  =  apparent  radius  of  sun  at  surface  of  the  earth  as  seen  by  satellite. 

As  outlined  in  the  Introduction,  the  different  J's  are  computed  for  a  given  V. 
In  the  case  of  ohmic  currents,  Jo^m  is  given  by: 


Johm<Xi> 


V(X.)  -  VlX-.j)  V(X.>  -  V(X.  +  1) 
R  +  R 


The  potential  V  is  then  varied  until  either  Eq.  (3)  or,  for  a  single  point,  Eq.  (2) 
is  satisfied. 


FORTRAN  LISTINGS 
Program  Units 

On  the  following  pages  are  listed  the  iterative  charging  program  SHADOW  and 
its  subprograms.  SHADOW,  itself,  is  only  a  blanket  program  directing  the  calls 
to  data  input  and  processing  subroutines  and  performing  data  output  directly.  The 
first  set  of  executable  cards  directs  the  input  of  plasma,  structure,  and  other 
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parameters,  with  the  option  of  calling  either  SBR  or  SPS  for  reading  structure 
data.  The  next  group  is  a  DO  loop  which,  for  a  series  of  structure  positions, 
calculates  voltages  by  calling  CNVG  and  then  prints  the  results. 

Subroutine  PLASMA  reads  the  values  of  number  density  and  temperature 
from  which  are  calculated  current  density  for  both  populations  of  each  species, 
protons  and  electrons,  and  reads  the  photoelectron  current  density. 

Subroutine  MODEL  determines  the  positions  of  the  structure  for  which  poten¬ 
tials  are  determined  and  stores  values  in  common  block  ZERO. 

Function  CNORM  determines  current  density  from  number  density  and  tem¬ 
perature. 

Subroutine  SBR  reads  the  dimensions,  conductivity,  and  number  of  nodes  of 
a  circular  structure  and  calculates  the  surface  area  of  each  node,  the  distance 
between  nodes,  and  the  conductance  between  nodes.  Solar  Power  Satellite  per¬ 
forms  the  equivalent  operations  for  a  square  structure. 

Function  CURV  determines  the  current  to  a  structure  as  a  function  of  potential 
and  position  on  the  first  node,  assuming  current  balance  to  all  but  the  last  node. 

Function  CURV  is  the  net  current  to  a  node  of  potential  VC  and  position  X. 

Subroutine  CNVG  finds  the  potential  on  each  node.  For  low  conductances, 
the  method  is  to  set  each  potential  at  the  zero  plasma  value  and  then  adjust  poten¬ 
tials  to  counter  current  flow,  eventually  reaching  equilibrium.  For  high  conduc¬ 
tances,  the  zero  value  of  CURV  is  found,  defining  the  values  of  array  V. 

Function  ZERO  finds  the  zero  value  of  function  FUNC. 

Function  DCDV  finds  the  derivative  of  CURV  with  respect  to  V. 

Function  CT  finds  the  total  plasma  current  density  to  a  node  of  potential  V 
and  position  X. 

Function  CJPE  is  the  photoelectron  current  density. 

Function  LNEX  determines  the  value  of  CHI(- V). 

Function  BS  determines  the  value  of  c. 

Function  SE  determines  the  value  of  a. 

Function  SI  determines  the  value  of  b. 

Function  FI  is  a  service  function  for  BS,  SE,  SI. 

Function  F  determines  the  fraction  of  sunlight  striking  the  spacecraft. 

Subroutine  AVG  is  a  service  function  for  function  ZERO. 


COMMON  BLOCKS 

For  high  efficiency,  COMMON  blocks  are  used  for  communication  between 
subprograms.  Thus,  it  is  essential  that  their  content  be  known  to  the  programmer. 
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Block 


Variable 


Purpose 


NP 

NP 

Number  of  nodes  (points) 

BOSTON 

SEP 

Separation  between  nodes  (km) 

X 

X  for  most  illuminated  node 
m 

III 

Service  index  variable  for  function  CURV 

V(112) 

Array  of  potentials  (volts) 

CON 

CON  (112) 

Array  of  conductances  between  nodes  (ohm-^) 

ALL 

Cl 

Ion  population  1  current  density  (A/m  ) 

CE 

Electron  population  1  current  density 

TI 

Ion  population  1  temperature  (eV) 

TE 

Electron  population  1  temperature 

CP 

Saturation  photoelectron  current  density 

TWOMAX 

CI2 

2 

Ion  population  2  current  density  (A/m  ) 

CE2 

Electron  population  2  current  density 

TI2 

Ion  population  2  temperature 

TE2 

Electron  population  2  temperature 

STRUCT 

AREA  (112) 

2 

Array  of  nodal  surface  areas  (m  ) 

GCON 

GCON 

Conductivity  of  structure  (ohm  1  m~*) 

MATER 

AA 

Material-dependent  secondary  and  backscatter 

BB 

emission  coefficient 

CC 

SHAPE 

SHAPE 

Value  is  1  for  spherical  probe;  0  for  plane  surface 

USE 

The  first  card  of  a  data  deck  specifies  the  distance  between  spacecraft  posi¬ 
tions  and  the  number  of  positions.  The  first  value  is  in  E  15.  6  format  and  the 
next  is  in  14  format;  the  distance  is  in  kilometers. 

Next  come  the  cards  specific  to  each  run:  the  cards  specifying  the  plasma 
conditions  and  the  structure  parameters.  The  first  card  is  a  title  card,  with 
alphanumeric  data  in  columns  1-57  and  the  spacecraft  shadow  potential  in  columns 
58-80  in  F22.  10  format.  In  the  listed  version  of  SHADOW,  the  shadow  potential 
is  ignored. 


The  next  card  lists  the  plasma  parameters:  ion  population  1  number  density, 

electron  population  1  number  density,  ion  population  2  number  density,  electron 

population  2  number  density,  and  photoelectron  saturation  current  density;  the 

3  o 

number  densities  are  in  n/cm  ,  the  current  density  is  in  A/m  .  Next  come 
the  four  population  temperatures  in  eV:  ion  population  1  temperature,  electron 
\  population  1  temperature,  and  so  on.  The  values  on  the  first  card  are  listed  in 

!  E15.  5  format;  on  the  second  card,  F15.  5  format. 

!  The  last  card  for  each  run  lists  the  thickness  (m)  and  conductivity  (ohms"' 

m  of  the  structure  in  El 5.  5  format,  and  the  width  or  diameter  (km)  and  number 
of  nodes  in  FI 5.  5  format. 


Sample  deck 


0E8UG=0UTPUTl 


- EKPEPHAL  3T - 

OIHENSIO'*  XL(112>  ,  F-LU12I.  VS*(112>,  VR(112>.  /  VSPI 1 121  , 

- 1  VWEtHttH - - - - 

COMMON  /*»P/N»  /80ST0N/SEP,  X,  III,  V(U2>  00029? 

- 1 - /ooHPeoxmei - vrllv9i>  cet ep - 

2  /ZERO'ACCUl,  ACCU2/VC/VC(112)  000010 

- 3 - Ff#WX23iei  9*2,  -m, -FEE - 000019 

4  'STRUCT/  AREA  J112I/5J1N/  SCOH 

- 5 - 'HATER/  W|  89  ,  ee - /3*R-»S/  SHAPE - 

OIHENSIOW  CUR  (1121,  CNT  (112),  TIT.E(S)  000020 

- B»fo-ap+g»evpf»ts  r— m— Ot  ■» - 

AA  =  1.3 

- 88-  --a. -a - 

CC  =  0.55 

8 - RE  18  8002.  BATA - 

CALL  MC35-  OiLTA,  N) 

- g  eewTi  wog - 

READ  (5,  ,0)  (TITLE (I),  1=1,6),  VAC 

- Ml - FORMAT  (  5-Altf  A-Er  P-22.101 - 

IF  (EOF(  5)  . AE.  0.  )  STOP  000920 

- WRITE  t*, All  -  <T-I-T-t"E<I»-|  -I-ti  6).  VAC - 

41  FORMAT  (LAI,  5A10,  AT,  FID. 21 

- e»tt-Pt-A-94A - 

SHAPE  =  ANi 

JOtteeT  ORTA  T-V  CA . StlLAT-;-  VOLT  A9E3  ~8W A  -SOW) RE  STRUCTURE - 

CALL  SP5 

COHHEMT-  IH  68T;  R- TO- CAtCULATE  THE-  VO-.TASES  OH  -A  -R-OUvO  STRUCTURE,  USE  STATEMENT 
CALL  SflR 

- V0  ■  ZERO  (Of, — •098,) - 

WRITE  (6,  701 

- 20-FORWA-T  -<tHOt-l-X-,— 5XefAtTfAEf-79)(^5A(t8€-AL  /  11X.  1SI1H-0,  24K, - 

1  65(1  M— )  /  3( ,  2HV3 ,  8X,  2rt(H,  5X ,  5HF(XH),  3X,  6HENERGT,  2X, 

- g  OHMET  SIR.,  6Xi  ?AXHr  5XT-5(F(XH)TTXr  lMVr-6XT— 44VVV0,-6«, - 

1  3HVVZ,  .X,  5HVVI/A,  SX,  6HSM  EPGT ,  2X,  8HMET  02.,  2X,  9H0AH.  CUR. 

- *~Vi - 


-O 


II  =  113  -  N’ 
08  1  I---t-T*- 


X  =  -135.2  ♦  -L3AT  (I-1)»3ELTA  -  0.5»SE» 

CALL-  0WV6 - 

XX  =  <  -  3F°*RL3AT(N»-1)'2. 

PF  »  «=•«  XXI - 

ENER5*  =  )  . 

8)  t  —  y  • 

DO  U  I  =  1,  H P 

IH‘  •  i -  - -  - 


X_  l  J  A  =  X  -  rL0AT (  J-l)  *SE* 

CJT  =  3 JT  »  WWh  <L(JM»APEA<» - 

CJ»!J)  =  CJR.  (V(JI,  <-  FL0AT{J-1)*SE») 

FF(.(J)  =  F(X.IJ))  - 

V’(J)  =  V<J('V0 

V5a(  Jl  -  V(l)  -  V(J»  — -  - 

Ir  (  V(l)  .ME.  V(HP|)  VV3*  ( J)  =  VSPU)  /(VUI-V(NP)) 

1=  t  J  i  £0,  HP) - 68  T8-  19 - 

rs((J)  =  <V(1)  -  V  (J*.  >  )  *C0N<  J) 

H  (99HE.H  .59,  9,1 - 59- T9  95  - - 


90-1*45 


Hits  PAGJJ  IS  B£SI  QUALITY  PSAaUCMIjS 

mo*  copy  masaisKKD  10  qdc _ 
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COMHON  /NP/N» 


2  /3ION'  GC3W 


100  F0»H#T  (-  15. 


/aosroM/SEP 


WHITE  (5,  1C11  M< »  JCON,  S,NP 


1  26MCON3u:TIi/ITr  MHOS/NETEHI  =  ,  £  12.  1  I 


2  r  1 0.  2»  13  X»  2  7HN'Ji9EP  OF  PJINTS  It)  MOHEL  =,  13  ') 


CON  <11  s  GCJN»T<<'S  /  SEP 


»*£*<!> 


ON<  1  J 


HiTUPN 


THIS  PAGE  IS  BEST  QUALITY  mCriOBL* 

.■^GUi»T-LV.aiSWa)TODDC 


- SUtHHHHIN:  P.«3H« - — 

COMMON  •  STRUCT/ARFl  <H2) 

- 1- - /Att/  oil  0 E,  TI,  T€t  i° - 

2  /TNOM1X/  CT?,  CE2,  TI2,  TE2 

- RS-At  KI,  <E - - - 

OUT*  KI»  <  E  /  1.24  F  *5,  5.31  E  »  6/ 

- reab  <9»  /■»>— oft— ore,  oci  o:g,  o»,  Tfr,  ?!?»  tc,  fee - 

70  FORMAT  (5  E  15.5  >  4  F  15.51 

- 61  ■  ONORN  (frl,  M,  X I ) - - - 

CI2  =  CNORM  (012,  1 12 1  KII 

- 0€  »  ONORM-tO-e,— T-~  >-«€) - 

CE2  =  CNORM  ( OE2 1  TE2,  <EJ 

- 0».  . . t  ,€ .  j - 

•MUTE  (6,  721  01,  IE,  012,  1E2,  ft,  TE,  TI2,  T->,  Cl,  CE,  CI2,  CF2 

- 1  ,  dp - 

72  FORMAT  (  30H3ION  NJMBER  OENSITY  (POP.  1)  =,  El?.  4,  5X,  35MELECTRON  N 

- 1  NUMBER  TENSITY  -(-POP,— 1M  «  -El-gi  4,  OX,  V - 

2  30H0  ION  NUMBER  TENSITY  <P0>.  2)  =,  E12.4,  5X, 

- 8*  '>5'MC'1E8T R BN  NUMBER  OENSITY  (  PBPt  -g)  -  ,  C  If  i  4  * - 

4  27M  01  ON  T  EM*  E  RAT  JR  E  (POP.  II  =,  F12.4,  5<, 

- 5-8*UCt€»T»«N-f eM»ERAMMRE---(P»»-,  t*  «|  E 12  1 4 ,  0X| - 4H<E(tt  / - 

6  27HOION  T  EM*  E  RATR JE  (POP.  ?»  =,  F12.4,  5«, 

- 7--31MC1E0PR0H- f£M»~(ATUR:  (POP,  g»  Fl2,4  / - 

8  31HOION  CURRENT  OENSITY  (POP.  II  =  ,  E12.4,  5X, 

- 5  35NEt€ttf  RON  TORRENT  8W1M  <PB»  I  t»  -  ,  0  2,41  OIH  ANfSP**) - / - 

A  31HGI0N  CURRENT  OENSITY  (F0>.  21  =,  E12.4,  5X, 

- 0  35MCtfcOIRON  CURR-ENT--BCNSITY  (POP.  2)  ■  ,  Cl  2.  4  < - 

C  32H0PH0r  0  EL:  C  TROM  CURRENT  IENSKY  *,  E12.  4  ) 

- RETURN - 

ENO 


- SUBROUTINE  NO  OCI  (OELTA,  N»  - - - 

COMMON  /ZERO/  ASCII,  ACCU2 

- READ  (5,  l«) - M6t.TR  ,  ~N - - - 

50  FORMAT  (  E  15.6,  14) 

- ACSU1  ■  8.  01 - - - 

ACCU2  =  C.01 

- MRITE  (>!■  5H  A 03 It  |  ACOUe,  OCLTR,  N - 

51  FORMAT  ( 1  Ml,  5MACCJ1=,  E13.6,  7H  ACCU2=,  E13.6,  OELTA=, 

- 1  -13,6,  >1  N»i  14  //» - 

RETURN 

- etrt— — — - — - 


IlHIS  PA®5  lS, 

pt  -.  .;i  COr  'i  ' 
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— rUMO-MON  ?!  <  *<  T)  — — - 

COHHON  7»A TE17  AA,  81!  CC 

- OINE-N9ION  A(3),  '3<!>,  C(3)-I~9-4-9» - 

OATA  A(LI,  A<2),  A  (  31  7  !M  3473F41  ,  .18440  08'  -  3 ,  -.5S9939S  7, 

-1 - 94t>,  9(2),  9(3)  7  -.1)  SSIS’E-S,  .1446707  -  3,  -,S715»19f-4  ft 

1  CC1I,  C(  2  ) »  “(31  7  •*6326013.  .9738416S-4,  .4775591E*1  7 

- IF  (4  ,»T,  8.) - 99-T9— t - 

OCH  =  FI;(V,A) 

— 0(g>->  FfKV,  )> - - — - 

0(3)  =  FI4V,3) 

- SI  »  rur,  0)  »  -so - - - - - - 

RETURN 

1  91  »  It  - — — - 

RETURN 
— gfrt - 


- F WOT  I0N--9-6— (-»t  'ft - 

CONNON  7  NATET  7  AA,  83.  CC 

- BIHfNSION-A45),  9f»» - 

DATA  A  (1 1  ,  A<2),  A(3)  7. -.8,24,  -.j2S6,  .3431  /, 

-1 - BUT  .  fl{  2  )  ,  9F3T7  .  663  y  -.-0-8frar  .l'r31-/ - 

IF  (V  .IT.  0.  )  G3  TO  2 

- IF  (T  . .  Ei  g)Ct)  39-TO-l - 

SE  =  FI(T,  B)  '  AA 

- ^ftufn - — - - - 

1  SE  =  FIIT.A)  •  AA 

— RETUPN - - — — - 

2  SE  s  a. 

- ,-^TURN  - 

EN1 


■I 

•i 


FUNCTION  "943  (  V  i  f-?1 - 

CONMON  /NATE 77  AA,  B8,  CC 

0INEN5I0N  A ( 3  )  ,  8433,  0(3),  0(3) - 

OATA  AUI|  A(2)  ,A(3)  7  - .  Ou  32) ,  15 1 27,  .  11 536/ 

-BATA  041)  ,  0(2)  ,9(3) - 7, 00019  ,.  0221  7, -.  000817 

OATA  Clll  ,C(’>  ,C«3|  7-. 02013  , -.  00901,.  2630467 

09  -»  --e,7  )H  «  130 - 

IF  <V  .L-.  0.)  RETURN 

if  e> .  > — *  »  go, - 

OU)  *  FTlUfjA) 

3(3)  *  Fjj v,3) 

■39  »  FI (  Ft1  |  ))  *—99 - 

IF  (OS  .  >T  .  1 . )  33  =  1. 

fB-  (09-tE-T  i  ).  ) - 39  «  0, - 

RETURN 
-frNB - 


%  .  .j  ',L"  * 

Y  _  A  G'O'E'  I 


ro  'Ji>u 


1 

I 
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I 

i 


— p:xc  hinot io<  lh:x  m> - 

COMMON  /SHXP :/  SH*»E 

If  w  .  >1  .  ***** - 60  T9-  -3 

IF  ( VT  .  .r .  -SOO.  I  60  TO  2 

— if-  tor  .lt  .  3. » — w-fe-1 - 

LNEX  •  1.  ♦  <r»SH»»E 

— r-:t  urn - — - 

1  ln-;x  *  e*»  c  vt  i 

Rr+UPN - - - 

2  LNEX  «  0. 

- RETURN - - — — - 

3  LNEX  «  VT  ♦  1.6*133 

- RETURN - — - — - - 

ENO 


pjNettwr-i-jp;  <vt*i - 

COMMON  /ALL/  Cl,  Cl,  TI,  Ti ,  CP 

OX  TX  OMRERE  *X  ./ - - - - - - — 

C  =  1. 


<SHX»g  .ET.  0»M-:pf> — C---3  .  o- 
if  cv  .;r.  o. >  ;o  to  i 

OJ»£  ■  0»|»F<<>»3 - 

RETURN 

-t-e^r  »-  op*-f«>a<i.  !»>«,»«»:  •  -e- 

RETURN 

- E<w - - - - - - - — - 


/SHXPE/  3 MX  PE 


i 


i 


- SUTROUTIw:  XVi  IV,  »i  0> 

o  -  i 

- V  ^-TV-»«TV  2. - 

RETURN 

- cNR -  - 


THIS  PACE!  IS  8i.5T  QUALITY  PEAGIICABIa* 
CO,-  i  i •  Ur-iilsKiS  TO  M>C 
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i 


REAL  LNE< 


common  i ci>  Tii  r : t  cp — /thomakf  CIg»  :!:e»  Tle»  T:e - 

cthvi  =  cjeivt  -  :jkv>  -  :i8S£tv»  -  cjsmv>  -  :jse<vj  -  :jpecv,x 

"H - 

C  JSE  (  V)  =  CF‘.NEX(V/TEJ*SE(/.  TE)  ♦  CE  2*LNE  X(  \l/r  ■_  >  )  *SE  (  V  ,TE  2 ) 

cjsitvt  *  et*t*e**-**n*+-9&  »>Tt'  ■  <  t-^»Me>*sn<»ria> - 

CJ8SE(VI  =  CE»LNEX(V/T£)  »8S  <V,T  E»  ♦  C  E2*L  NEXT  VX  r  EE  )*  BS  TV  >  T  E2  ) 

CJ  I<¥>  -  3I»-N£X<-T.>TI>  ♦  Cl  Xt-V/TI2> - 

C  J€  f  V  >  *  CE’.NEXTV/TE)  ♦  CE  ’  M.NET  <  V/TE  2) 

IF  <A8ST»-0.j> — .ST.  O.S>  ft o  ro-g - 

CT  =  (CTK1.I  -  CTIU.M’V  ►  CTi:0.) 

-RrTNRN - 

2  CT  =  CT1 J  V  ) 

- RETURN - 

ENO 


- FUNOTIOV  DCOV 

COIMON  /CD N/ 

<V'T  II - 

CONUO) 

/STRUIT/  AREA  (t0> 

f  N  P/  NP 

CC  =  COKL ) 

1  IF  a  .n-:.  N>l 

GO  TO  2 

GO  TO  3 

- 2  08  =  CONTI >  *  COMT I-l> - 

3  OC J V  *  5i  ♦  »REAtI»*(CT(»*:,  01,0.  I-CT(V,0.>)*1.'3. 

- R£-HN»N -  - 

ENT 


FUNCTION  TER)  TFJN%  X) - 

COMMON  TT;R0'«CCU1,  ACCJ2 

<t-  -  f. - — — - - 


VI  =  0. 

— r?  -  -etmc.  — - - - - - — — — — — — — - 

6  c  =  fun:  <  v,  x  i 

- IF  <T*i»<Vl-m — rfrr.-- AOOU1)  .AND. — HWW  rtEi  AICUeH  BO  TO  3 

IF  (C  .. 

CO  TO  b 

j  A1  •  »  A 

o CP.  L  wTJ 

GO  TO  ft 

s  f-"»o  «  v 


u.  t  GO  TO  A 
♦  ¥r  v ?!  VET 


**1  Vl«  »»> 


RETURN 


i 


l 
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r  ■ 

•k- 


,  ,  jAUU’i  d 


SUlROUTt  HE  00*6 - 

extepnal  :t,  c-Ml,  cup# 


/60STOH/S4Pt  X,  II, 


/2ER0/ACCU1,  ACCU2/YC/ VCtlOl 


/strict/  »re»(io)  'con /  comioi 
/tr  oh  ax  '-<?*-?»  e€g>  TteiTee - 


IF  (C  ONI  1 1  3.1  CO  TO  12 


P  =  1  ./CON  (II 
00  It  I  «■  if  HP- 


IF  Cjlll-ll  .FI.  8.1  GO  TO  12 


- ft-.— g— 1  ,  ,N)>N-<-F-  t  | - 

11  CONTINUE 

- ¥<11  ■  Z€*(H-Mr*H - 

IF((CT(/<1  »*j.  01,0.1  -CT(#(1>  ,  0 . 1 1  •ARFA(1)*10  0.*T  .LT.  ».0>  63  TO 

- »— tt - 

12  VS1  =  l.E»  103 

00  1  I  » 


l  V  (  1 1  =  ZE«  (CT,  X-FL0AT(I-1I*SEP1 

>  continue - 


NN  *  NP 

m  -3 — t— » — tr  +*— 

WC(I)  =  V'll 
- H— — t- 


»di»  =  yiin  -  curl i v i it  i .  x-FLOAT(it-i)»sr»»/oco#(  « in,  in 

»  o.j - 

II  =  NP.l-I 

*CtI{»  «T(-HI - 

V  ( I II  =  V(  T  T I  -  CURL  ( tf(  1 1 1  i  X-FL0AT(II-l»*SE»»/0C0iM#(in,  III 
*"0;3 - 


3  CONTINUE 

- - 

00  4  1*1,  NP 

— * - #0  *  #q **m — vc<n>»-»g 

IF  (VS  ._(.  tCCU2*»2l  GO  TO  E> 

— - IP«»9  tSfi  Pit  AOOJei  GO  ID  »- 

00  6  I*  1,  N3 

- 9 - #<I>  «  *Oi  I  ) - 

VS1=VS 

- GO  TO  2 - 

5  10  5  I  =  1,  IP 

— ‘-5 - #111  ■ — WH+-I  (OIDini - 


VSl=VS 

- GO  TO  t - 

6  CONTINUE 

- <?c  rw»»f  -1 - 

10  CONTINUE 

- #U>  •  >:1011UIWt  <1 

PETU®N 


-  -7 


.-i 
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i 


COMMON  fNP/N*  /CON' G  ONU  0) 

i — - f-tostyw/'h >f  -to  i  19' - 

Z  /STRUGT/  AREAI112) 

IF  t!  sHr.  t*  SO  f)  1 - - 

CURL  =  CT(W,<)»ARrAC)  ♦  £¥-0C<2)  )>C0N(1> 


v ( r  »  =  U<  I-l  )-IMr-2)>»C0N<t- 2>/C0N(I-l»  f  V  ( I-  1 »  *  CT<  If  <  1-1 )  » 


1  CONTINUE 


M  I  I  ‘  li  IF 


RETUPN 


COMMON  /**>/  IP  /BOSTON/  SE® 


onf  g:on 


OINTE G«,J  )  :  AINTEIWY)  -  a:NTEG(X> 


100  FORMAT  <’  E  LE.5,  ’  r  1=>.5> 


WRITE  (if-  10U  T  M<  f  GCON,  R,  NR 


■ARE  A  <  NPI  -  AREAUT 


RETURN 


*T.  •  '  - 
T'  .  •  A  S 


l  '  Ijt 


-  •  ;/cat:  rxJ&nfiiX* 
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